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Abstract. The Particle Number Projected Generator Coordinate Method is formulated for the pairing 
Hamiltonian in a detailed way in the projection after variation and the variation after projection methods. 
The dependence of the wave functions on the generator coordinate is analyzed performing numerical 
applications for the most relevant collective coordinates. The calculations reproduce the exact solution 
in the weak, crossover and strong pairing regimes. The physical insight of the Ansatz and its numerical 
simplicity make this theory an excellent tool to study pairing correlations in complex situations and/or 
involved Hamiltonians. 

PACS. 74.20.Fg BCS theory and its development 
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1 Introduction 

The measurements of Black, Ralph and Tinkham of 
discrete level spectra and spectroscopic gaps in nanometer 
Al isolated grains were interpreted as evidence of the su- 
perconductivity phenomenon. To understand the physics 
of such ultrasmall grains a great deal of theoretical effort 
was devoted to study such systems starting from grand 
canonical (BCS) and canonical ensembles as well as 
very sophisticated theories j4j of the pairing Hamiltonian, 
• see ref. IB] for a review. Later on the exact solution J3, 
of this naive model Hamiltonian was rediscovered. Some 
, others studies treat the aspect of thermodynamic proper- 
' ties |7||H| while others, see for example the question of 
' persistence of pairing correlations above the BCS critical 
i temperature is addressed. More recently some analytical 
' results in special regimes have been obtained ^Uj. The 
, main issue of all these studies is the proper description 
■ of the crossover between the few electron regime and the 
I bulk one. To analyse this crossover several properties can 
. be computed as a function of the mean electronic level 
' spacing d (or the number of electrons N) that character- 
izes the transition from one regime to the other. One of the 
findings of these studies was that the strong phase tran- 
sition predicted in a grand canonical study was absent in 
more advanced theories as well as in the exact solution. 
In the BCS approach superconductivity is not possible for 
all d (N) breaking down at a critical d value. This break 
down is number parity dependent and indicates that quan- 
tum fluctuations are not treated adequately by the BCS 
wave function. The knowledge of the exact solution for the 
simple-minded pairing Hamiltonian does not diminish im- 
portance to the theoretical approximations developed for 
the study of that Hamiltonian, see ^Tj for a review. These 



approximations are very general and allow the study of 
more sophisticated Hamiltonians for which no exact solu- 
tion exists. The use of exactly solvable Hamiltonians|12|. 
on the other hand, is very practical since it allows to check 
the accuracy of different approximations in the limiting 
cases represented by those Hamiltonians. 

In a recent paper U3| we have proposed a new ap- 
proach to study superconductivity in finite systems, namely 
the Generator Coordinate Method (GCM) based on 
particle number projected BCS wave functions generated 
in a suitable way. In that paper the GCM approach was 
applied to superconducting grains described by the Pair- 
ing Hamiltonian and it was shown to provide an accurate 
description of these systems in perfect agreement with 
the exact Richardson solution. The purpose of this pa- 
per is two-fold, first, to present a detailed derivation of 
the relevant formula as well as the way to solve the Hill- 
Wheeler (HW) equations and, second, to analyse different 
generator coordinates in the context of pairing correla- 
tions. The derivation presented is comprehensive enough 
to allow for the application of the formalism to other pair- 
ing Hamiltonians. Furthermore since our theory is very 
general and not constrained by any requirement can be 
applied to more complex systems. As a matter of fact we 
have performed preliminary studies with the most general 
pairing Hamiltonians proposed in |12) and the results |15| 
are of the same quality as the ones presented in this in- 
vestigation. Finite temperature effects are not considered 
in the present study. 

In sect. 12 we derive the general formula of the GCM. 
In sect. Owe discuss the different coordinates to be used 
in the calculations. The convergence and other issues con- 
cerning the numerical solution of the HW equations is 
analysed in sect. 01 Finally in sect. El the whole formal- 
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ism is applied to study superconducting grains. The pa- 
per ends with the Conclusions and some numerical aspects 
discussed in the Appendices A and B. 



2 Theory 

The pairing Hamiltonian used in most calculations is given 

by 

N N 

H = ekcl^ Ck,u 4+4~Cfc'-Cfe'+ (1) 

k=l,u=± k,k'=l 

where k+ (fc— ) labels the single particle level (time re- 
versed) with energies and Cfc , cj, destroys and creates 
electrons in their respective states. The interaction con- 
stant G is taken as Ad with d the level spacing and A 
the BCS coupling constant whose value for Al is 0.224. 
The single particle energies for simplicity take the val- 
ues £k = kd. The number N of electrons is equal to the 
number of levels and in the ground state they form N/2 
Cooper pairs, so one works at half filling. This Hamilto- 
nian allows the discussion of the crossover between the 
strong-coupling regime {d/A <C 1) that represents large 
grains and the weak-coupling regime {d/A ^ 1) for small 
grains, in terms of the quantity d/A — 2sinh(l/A)/iV with 
A the bulk gap, or equivalently in terms of the number of 
electrons N. 

The simplest way to deal with pairing correlations is 
provided by the BCS theory ^Q. Its Ansatz is given by 
the mean field wave function 



N 



Uk ■ 



(2) 



fc>0 



The variational parameters Vk are related to the proba- 
bility to find two electrons in the level k. The parameters 



Uk are given by 



1. The spontaneous particle 



number symmetry breaking mechanism implicit in Eq. El 
enlarges the available variational Hilbert space making the 
BCS approximation, in the case of large particle numbers, 
a very good one. The BCS state (jSJ, on the other hand, 
undergoes strong particle number fluctuations and for fi- 
nite systems like metal grains, the Ansatz (jSJ is unreliable 
and misses essential features. To correct this failure it is 
necessary to develop the BCS formalism in a canonical 
ensemble, where the particle number is fixed, rather than 
in a grand-canonical one. The restoration of the particle 
number in the BCS context was introduced by Dietrich 
and Mang |17j in a nuclear structure context and it was 
applied for the first time to superconducting grains by J. 
von Delft and F. Braun |18II19| . The projection method is 
based on the Anderson formulation of superconductivity 
i20 where projection onto good particle number is pre- 
sented as an integration in the gauge variable (/?, 



\BCS) 



N 



2ti 



-iip/2 



Uk + e 



(3) 



We assume the number of particles N to be even, the odd 
case is considered in appendix lA.ll The formulation of the 
particle number projection (PNP) can be done in several 
wavspi]. Very compact formula are obtained in terms of 
the residuum integrals 1171 defined by 



1 f^"" 
27r Jo 

N 



'i(M-2m)ip/2 



X n {^-'^"ul + e'-'\l) (4) 

This definition holds for indices ji , ■ • • , ju such that jk ^ 
jp for all k and p. The integer M is simply a counter of 
the j's involved. In case that two or more indices are equal 
we define the corresponding residuum integral as zero. All 
expectation values can be easily calculated in terms of the 
residuum integrals. As an example we evaluate the ma- 
trix element m{BCS\BCS) n, direct substitution of Eq.O 
provides 



{BCS\BCS)^ = / ^ \{{e-^^"ul + e^-l\l) ^ i?". 

(5) 

In the same way the projected energy is given by 

N {BCS\H\BCS) N 



N 



E 



N 



n{BCS\BCS)n 



= 2E - f ) § -G±u,v,UkVk?^. (6) 



The PNP energy, as the BCS one, depends only on the 
variational parameters Uk,Vk- Minimization with respect 
to these parameters leads to a set of N coupled non-linear 
equations 

2{ek + Ak)ukVk - Ak{ul - vl) = 0. (7) 
The quantities efe, Ak and Ak are defined by 

4 = (6, - G/2) |i , Ak^GY C (8) 



A, 



N 

E 

3 



R 



kj 



r\' R{ Rt 



Rq 



N 



R. 



- Rf 



kji 



tdJ^ Tjk T)k 
ri^ — -Kg 



(9) 



The set of equations ^ resembles the ordinary BCS 
equations. In that equations, Ak = Q and = e — Gvf. — ^. 
The Lagrange multiplier /i takes care, on the average, of 
the particle number conservation. Notice that in the pro- 
jected equations the fields Ak appear in addition. The so- 
lution of Eqs.Q defines \BCS)n- In the literature jlH] 
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this is usually called Projected BCS (PBCS) theory. De- 
tails of how the set of Eqs.(|2j) is numerically solved are 
given in appendix IA.2I 

To include additional correlations we consider a gen- 
eral superposition of different projected-BCS wave func- 
tions, 

l-^w) - y rfC fiO \bcs{0)n 

X n(e''^^'"fc(0 +e"^/'«fc(04+4J|->- (10) 
k 

The new wave function is based on the Gener- 
ator Coordinate Method (GCM) developed by Hill and 
Wheeler in Nuclear Physics It has been also used 
by Peierls, Yoccoz and Thouless among others to 

deal with the restoration of symmetries in mean field ap- 
proaches as well as to deal with a variational approach to 
collective motion. It has also provided a variational deriva- 
tion of the Random Phase Approximation ^24,. The coor- 
dinate ^ refers to any parameter on which the BCS states 
may depend parametrically. In this way the superposition 
state Itf^Tv) takes care of the fluctuations associated to the 
parameter ^. In principle, the variational quantities are the 
weights /(^) and the occupancies Uk{£,) , Vk{£,) and should 
be determined invoking the variational principle. The final 
equations, however, result in an integro-differential set of 
equations very complicated to solve. In consequence some 
assumptions about occupancies are needed in order to fa- 
cilitate the numerical implementation. If we assume that 
the quantities Uk{£,) , Vk{(,) are known (see below) one deals 
only with the problem of calculating the weights /{(,)■ This 
is accomplished by the Hill- Wheeler (HW) equation 

J d({n^^> -m^^>)f{0 = o. (11) 

Ti is the Hamiltonian overlap, defined by 
H^^,=n{BCS{0\H\BCS{^'))n 

3 

-G J2 M^')^^iOMOv3i^')RH^^^') (12) 

and Af the norm overlap 

U^^, ^N{BCS{0\BCS{e))N=R°o{i,^'). (13) 

This equation is very similar to Eq. Q . The residuum 
integrals have now been generalized by 

X [] {e-''''^Uk{S,)uu{i) + e'^'\u{Ovk{i)) (14) 



In App. rSl we discuss the way to calculate these integrals 
and the Hamiltonian overlap of Eq. H12|) . It is important 
to notice that the solution of the HW equation provides 
not only the ground state but also the low-lying collective 
states. 

As we mentioned above, in the HW equation the quan- 
tities Uk{^),Vk{^) are supposed to be determined before- 
hand. They are usually fixed by the way the projected 
wave function \BC'S{S,)) n is calculated, namely, whether 
\BCS{^))n is determined by projection after variation 
(PAV) or variation after projection (VAP). In the former 
(PAV), the occupancies are determined by the symmetry- 
violating wave function \BCS), i.e., by solving the ordi- 
nary BCS equations. In the VAP case the occupancies are 
given by the solution of the variational equations Eq.l(7I). 
In the BCS framework the VAP approach is known as 
PBCS. Obviously the VAP method is more involved but 
it is a fully self-consistent method that provides better re- 
sults. We shall denote the first method GCMPAV and the 
second one GCMVAP. 



3 SELECTION OF THE GENERATOR 
COORDINATE 

The generator coordinate ^ is quite general and its selec- 
tion is motivated by the physical problem. The BCS wave 
functions depend parametrically on the generator coor- 
dinate, its selection is therefore strongly related to the 
ways we have to characterize the wave function. Though 
there are many ways to choose the generator coordinate, 
we think that for the BCS case there are three relevant 
ones: The gap parameter A, the Lagrange parameter fi as- 
sociated with the particle number of the BCS wave func- 
tion and AN"^ = {BCS\N^\BCS) - {BCS\N\BCS)^ , the 
fluctuations on the number of particles of the BCS wave 
function. 

Instead of using directly the gap parameter as a co- 
ordinate it is numerically easier to generate BCS w.f.'s 
with different gap parameters by solving the correspond- 
ing BCS (PBCS) equations for different values Gtriai of 
the strength constant G. This method is easy to imple- 
ment and very efficient. The second method is the sim- 
plest one. Now the generator coordinate is the chemical 
potential /z which in the ordinary BCS equations is used 
as Lagrange multiplier to fix the mean value of the particle 
number in the grand-canonical ensemble. In our case we 
solve the BCS equations for fixed fi and the use of different 
fi values allows to generate wave functions \BCS{fj,)) with 
different average particle number. The fact that \BCS{^)) 
does not have on the average the right particle number 
does not matter since later on we project on the right 
particle number. In this case one is looking for the fluctu- 
ations in the position of the Fermi level. 

The last method, finally, considers fiuctuations around 
the uncertainty in the particle number AN"^. Now it is 
necessary to add a constraint to fix a given value of AN'^ . 
This is done by using the modified Hamiltonian H' = 
H — jiN ~ ii2AN'^ , where the parameter ^2 guaranties that 
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the constraint is fulfilled. In principle we have set up six 
variational methods (PAV and VAP versions of each coor- 
dinate) but only five are feasible, because the VAP version 
of fj, (by construction) is not possible. Although the ulti- 
mate test of the quality of the selection of the generator 
coordinate will be the eigenstates of the HW equation it is 
interesting to have a look on the Hilbert space generated 
by the different coordinates. The diagonal elements of the 
matrix H^^/JV^^ of the HW matrix, Eq. are the pro- 
jected total energies En{^). This quantity is related to 
the condensation energy (CE) by EconiO = En{£.) — Ep 
with Ep the uncorrelated energy of the Fermi sea, i.e., 
Ep = — GN/2. Though, as mentioned above, we 

have five ways to generate w.f. depending parametrically 
on ^ we shall concentrate in this section on the three PAV 
cases corresponding to the three different coordinates un- 
der study. In Fig we display EconiS.) as a function of 
the corresponding generator coordinate and for different 
particle (level) numbers to cover the full range from weak 
to strong pairing regimes. For simplicity we plot only the 
curves for grains with an even number of particles. Let 
us first discuss the coordinate G trial- It is obvious that, 
for each number of particles TV, a critical value Gc{N) of 
Gtriai exists such that no superconducting solution of the 
system is found below it. In panel a) we show the CE 
versus Gtriai — Gc{N). We find a parabolic behavior with 
the vertex moving to larger values of Gtriai — GdN) as 
the particle number decreases (as one would expect). The 
curves get softer with decreasing particle number with the 
curve N — 20 being specially soft. We also find that the 
value of the CE in the minimum is larger (in absolute 
value) as the particle number increases (as one also would 
expect). 

In panel b) the quantity EconifJ') is plotted against 
M ~ fJ'BCsiN), Hbcs{N) being the chemical potential of 
the BCS equation for the corresponding case. Because 
of the particle-hole symmetry of the model the subtrac- 
tion of fJ-BCsiN) provides symmetric curves around /i — 
A*scs(-^) = ^. For large N the solutions are parabola 
like curves which soften with decreasing particle number. 
For > 40 we find superconducting solutions for all /i 
values. This is not the case for N = 20 where for certain /i 
intervals we do not obtain any solution for the BCS equa- 
tion, see below for more details. This is not surprising be- 
cause the standard selfconsistent BCS equation does not 
provide a correlated solution in this case, see below. Lastly 
in panel c) Econi^N'^) is plotted against AN'^. Here we 
also obtain a parabolic behavior similar to the case a) with 
the difference that the minima shifted to large AN'^ cor- 
respond to the large particle numbers. The CE gets softer 
with larger particle number as one would expect. 



^ It is important to notice that in this case the num- 
bers 20, 40, etc correspond only to the number of levels 
and not to the number of particles of \BCS{^)). Since we 
work without constraint on the particle number in general 
{BCS{^J,)\N\BCS{^)) / iV. On the other hand since we 
are projecting on the particle number the wave functions 
\BCS{ij.))n correspond to a system with N particles. 
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Fig. 1. Projected condensation energies, in units of the bulk 
gap, as functions of the different generator coordinates in the 
PAV approach. 



It is clear that the energy minima of the different co- 
ordinates provide an approximation to an unconstrained 
VAP calculation. In Table^we have summarized the min- 
ima of the parabola as well as the VAP values and the ex- 
act ones. We find that all three coordinates do a good job 
for large particle numbers and that big differences appear 
for small particle numbers, i.e., in the weakly correlated 
regime. We find that in general and at this level the co- 
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Table 1. Condensation energies, in units of A, predicted by 
PAV, VAP and exact calculations. 



N 


20 


40 


86 


172 


400 


Gtrial 


-1.7716 


-1.8194 


-1.9053 


-2.3566 


-3.4192 




-0.7864 


-0.9272 


-1.4925 


-2.0392 


-3.4625 


AN^ 


-1.1438 


-1.3654 


-1.6906 


-2.2227 


-3.5564 


VAP 


-2.0625 


-2.2441 


-2.4015 


-2.5428 


-3.6551 


exact 


-2.2026 


-2.5284 


-2.9403 


-3.5322 


-4.8891 



ordinate Gtnai is the most effective followed by AN"^ and 
11. Of course this does not mean very much since the con- 
figuration mixing calculations will change these results. 

Let's now analyse the wave functions generated with 
the different coordinates. To a given value of a coordi- 
nate, let say ^o, corresponds a wave function \BCS(S,o >■ 
A simple way to characterize the physical content of this 
wave function is by the associated gap parameter A(^o) = 
^J2k ^k{^o)vk{£,o)- In Fig. 121 we have represented the gap 
parameter A{^) associated to each wave function as a 
function of the coordinate ^ used to generate it. In panel 
a) we show the results for Gtriai ■ Of course the G entering 
into A is the one of the original Hamiltonian, see Eq. , 
independently of the Gtriai used in the calculations. Tak- 
ing into account the expression of A we expect, in first 
order, a linear behavior with Gtriai and this is what we 
obtain. In general a very broad range of gap parameters 
is covered, which is the reason why the coordinate Gtriai 
can be considered equivalent to the gap parameter A. The 
case of the coordinate fi is considered in panel b), where 
we represent the corresponding gap parameter as a func- 
tion of n—^gfjg. We find an oscillating behavior of A with 
/i due to the symmetry of the model. Notice that the scale 
of the y-axis depends on the particle number considered, 
see the figure caption. For fj, = kd, i.e., at the single par- 
ticle energies e^, we find maxima and for ji = k{d + 1/2) 
minima. The period and amplitude of the oscillations de- 
crease with growing particle number because in this model 
d ~ 1/N. For > 40 we obtain superconducting solutions 
for all /X values, in particular for /i — ^ibcs, i-e., for the 
selfconsistent BCS equation. For N — 20, however, we 
observe that at and around fi = k{d + 1/2) we do not ob- 
tain correlated wave functions. As mentioned above this 
behavior is in agreement with the fact that the selfcon- 
sistent BCS solution does not have correlated solutions 
in this case. The situation is further illustrated in Fig. O 
for the iV = 20 case. In the weak pairing regime we only 
find solutions for /i values corresponding to the hatched 
regions around a given level. In the region around the 
level fc, the number of particles of the BCS w.f., i.e., the 
expectation value {BC S{ii)\N\BC S{ij)) , varies in a con- 
tinuous way from 2(k — 1) to 2k. For example for the case 
of iV = 20, i.e. k = 10, the BCS w.f. around the level 10 
have average numbers of particles ranging from 18 to 20. 
In general from these w.f.'s it is always possible to project 
to 20 particles. In the regions between the hatched regions 
no BCS solution is found but only the Hartree-Fock (HF) 
one. The numbers of particles are obviously integer num- 
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Fig. 2. Dependence of the order parameter A on the generator 
coordinates Gtriai , fJ- and AN^ . In panel (b) the y-axis scale 
applies only for N — 20, for A*' = 40 the y-axis covers the 
interval 0.6 - 1.4 and for N = 86, 172 and 400 the interval 
0.9 - 1.1. 
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14 . 



8 



d 



Fig. 3. Sketch of the regions of weak and strong pairing for 
A'^ = 20. The numbers on the right hand side correspond to 
the labels of the levels while the ones on the left hand side to 
the average number of particles of the BCS w.f. at the corre- 
sponding II. 



bers, in the example displayed these integers are 16, 18, 20 
and 22. To project to 20 particles from these HF w.f. is 
only possible for 20, in the other cases the w.f. is zero. 
This fact explains the curve corresponding to = 20 in 
Fig. n From all regions where no BCS solution is found 
only the one between fc = 10 and fc = 11, corresponding 
to a HF solution with 20 electrons with zero condensa- 
tion energy, survives. In Fig. ^ this region is represented 
by the straight hue around n = hbcs- For TV > 40 this 
is not the case and we always find BCS solutions. The 
hatched regions of Fig. |21 correspond in this case to strong 
correlations and the white regions to weak ones. 

Finally, in panel c) the gap parameters corresponding 
to the AN"^ generator coordinate are plotted. The behav- 
ior is again linear, as for the coordinate Gtriai, but the 
range of the gap parameters involved in each wave func- 
tion is the opposite one. In this case we obtain for small 
particle numbers a much larger range than for large par- 
ticle numbers. 



4 NUMERICAL SOLUTION OF THE 
HILL-WHEELER EQUATION 

For our purposes solving the HW equation, Eq. ^2 is 
equivalent to the diagonalization of the Hamiltonian in the 
nonorthogonal basis of the generator states \BCS{^)) n ■ 



The usual procedure to deal with this equation TT in- 
volves two diagonalizations. In a first step, the norm over- 
lap A/'^j' is diagonalized 



(15) 



with the functions Mfe(C) forming a complete orthonormal 
set in the space of the weights f{^). Its eigenvalues are 
never negative, Uk > 0, because the matrix J\f is definite 
positive. We shall keep the Uk{^) with nonzero eigenvalues 
corresponding to the linearly independent states. In prac- 
tice and due to numerical reasons we restrict the Ufc(^) to 
those with eigenvalues larger than a tolerance e. For each 
of these functions there exist states |fc). 



|fc) 



1 



/Uk 



dCukmBCS{0)N, 



(16) 



called the natural states, which span a collective subspace 
He- In a second step the Hamiltonian H is diagonahzed 
in this space 



J2{k\H\k')gk' = Egk 



(17) 



with 



{k\H\k') = 



1 



dCdC'u*(OWcC'"fc'(C')- (18) 



^/nkTik' 

The HW equations provide a set of wave functions. 



If 



Nl 



E 9l\k), 



(19) 



and energies Ecr labeled by the index cr, the lowest one cor- 
responding to the ground state and the others to excited 
states. In this work we are only interested in the ground 
state. Taking into account Eg. 1101 and Eg. 1191 one obtains 



9k 



UkiO- 



(20) 



Since the wave functions \BCS{^)) n are not orthogonal, 
the weights /(^) cannot be interpreted as the probability 
amplitude to find the state \BCS{(,))n in I^'at). It can be 
shown, however, that the functions 



Gio= E 9kUkio 



(21) 



are orthogonal and that they can be interpreted as prob- 
ability amplitudes. 

For numerical purposes all the expressions above in- 
volving integrals have to be replaced by sums discretizing 
the space ^. In this form one deals with matrix equations 
easier to handle. The question that immediately arises is 
how to determine the optimal ^-mesh to be used in the 
calculation. The border values of ^ are determined by en- 
ergy arguments since the probability of mixing high-lying 
states is very small. The ^-coordinate intervals used in the 
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calculations are given in Table El The calculations depend 
furthermore on the mesh step used in the discretization. 
This parameter is chosen as to optimize the calculations, 
i.e., we take the largest mesh that includes all states with 
relevant information. This parameter is also related to the 
required accuracy. In the calculations performed we have 
not attempted to reproduce the exact results up to an un- 
usual accuracy. In Fig. ^ the convergence of the conden- 
sation energy, as a function of the number of mesh points 
used in the calculations, is shown for different numbers 
of particles and for the coordinates G, fj, and AN'^. We 
observe that the number of mesh points needed for con- 
vergence depends on the ^-coordinate and on the number 
of particles. The coordinate Gtriai, see top panel, provides 
the best convergence of the three calculations. For large 
particle numbers a very good convergence for relatively 
few mesh points is found. For small numbers of particle 
one has to go to larger mesh points to find the plateau. For 
the fj, coordinate the situation is the reverse one, i.e., one 
finds earlier convergence for small numbers of particles. 
Finally, the situation for AN'^ is something in between 
the two former cases. We find that to reach convergence 
in energy 80 mesh points are sufficient for all coordinates. 
This is the number which we will use in all following nu- 
merical applications. 



Table 2. Initial and final values of the generator coordinates 
used in the calculations. 



N 


20 


40 


86 


172 


400 


G.(meV) 


0.44 


0.18 


0.07 


0.03 


0.01 


GfimeV) 


0.80 


0.30 


0.17 


0.06 


0.04 


Hi{meV) 


-3.00 


-2.00 


-2.00 


-1.00 


-1.00 


fif{meV) 


3.00 


2.00 


2.00 


1.00 


1.00 


ANj 


0.00 


0.00 


0.00 


0.00 


0.00 


AN} 


8.00 


12.00 


16.00 


16.00 


32.00 



A further check concerning the convergence is the num- 
ber of natural states kept in the calculations. Since many 
of these states are linearly dependent some natural states 
\k) will have a vanishing norm and must be excluded. 
In the calculation only those natural states with a norm 
larger than a given tolerance e are kept. For a given tol- 
erance we take as many states \k) as needed to reach a 
good plateau. Now we analyze the energy convergence as 
a function of the number of natural states kept in the diag- 
onalization of the HW equation or equivalently of the tol- 
erance of the calculations. This is shown in Fig. [3 for the 
three coordinates and for grains with different numbers 
of particles. In the fi and AN'^ coordinates we find that 
for tolerances smaller than lO^^*' linear dependent states 
are introduced in the calculations providing unrealistic en- 
ergy values. The interesting point is the nice plateau found 
for larger tolerances. The tolerance of 10~^° corresponds, 
typically, to around 15 linearly independent states. The 
coordinate G is in this respect somewhat different. One 
observes that for tolerances of up to 10^^^ one still gets 
linearly independent states, which obviously correspond 
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Fig. 4. The condensation energy in units of A for the three 
coordinates as a function of the number of mesh points used 
in the calculations. The energy scales correspond to the N=20 
case, the other curves have been shifted in order to make the 
figure readable. The shifts are 0.28, 0.58, 1.02 and 2.48 for 40, 
86, 172 and 400 particles respectively. 



to highly excited states that do not affect the energy of 
the ground state. This tolerance typically amounts to 20 
linearly independent states. From this respect we conclude 
that if one is interested in excited states the coordinate G 
is more effective than the other ones. 

The diagonalization of the HW equation in the case of 
the /i coordinate requires some comments. As mentioned 
above in the weak pairing regime, for = 20 for exam- 
ple, and for several /i intervals one does not find supercon- 
ducting solutions. This circumstance, as explained above, 
shows up as "missing points" in the CE curves. These dis- 
continuities do not affect however the solution of the HW 
equations, since these points have norm zero and do no 
mix with the other states. 



5 APPLICATION TO SUPERCONDUCTING 
GRAINS 

In this section we present a systematic study of properties 
of superconducting grains. The results of the GCM calcu- 
lation of different quantities are compared with the PBCS 
approximation and the exact Richardson solution. 
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Fig. 5. Same as Fig0]but as function of the tolerance e. 

5.1 Ground State Condensation energies 

Condensation energies characterize the presence of pairing 
correlations. The crossover between superconducting and 
fluctuation dominated regimes can be described through 
this quantity. 

As in the former cases the condensation energy Econ is 
defined as the difference between the total energy in the 
corresponding approximation and the energy of the uncor- 
related Fermi sea. For example in the GCM approaches it 
is given by Econ = Ecr=o — Ep, see Eq. ^1 and below. 
This quantity is displayed in Fig. for even grains (up 
to 600 electrons) and odd grains (up to 601 electrons) as 
a function of the particle number N. In both plots we 
give numerical results for the approximations discussed 
above, BCS, PBCS and the GCMPAV and GCMVAP ap- 
proaches. The GCMPAV results are presented for the co- 
ordinates G trial , M and AN"^ and the GCMVAP for G trial 
and Z\iV^. The grand-canonical (BCS) calculation of Econ 
predicts vanishing correlations in the few-electron regime 
in the even and odd systems. The PBCS condensation 
energies, on the other hand, though always negative pre- 
dict an unrealistic sharp crossover between the fluctua- 
tion dominated regime and the bulk which is more pro- 
nounced in odd grains. This artifact is not present nei- 
ther in the GCM approaches nor in the exact calcula- 
tions. The simpler GCMPAV approaches already predict 
a smooth crossover for odd and even grains. The more in- 
volved GCMVAP approaches not only predict a smooth 
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Fig. 6. Condensation energies versus the number of particles in 
different approximations and the exact results. Upper (lower) 
panel for even (odd) systems. 



crossover but their predictions coincide with the exact re- 
sults. Concerning the GCMPAV calculations we find that 
the fi coordinate is the most effective of all of them fol- 
lowed by the AN^ one. Paradoxically the calculation with 
the ^ coordinate is the simplest one from the numerical 
point of view. 

The reason why the /i coordinate is the most successful 
one is probably due to the fact that using this coordinate 
one has the right inertia parameter for the rotations in the 
gauge space associated with the operator N. As a matter 
of fact this was demonstrated by Peierls and Thouless ^H] 
in the context of the translational invariance and a PAV 
approach by the double projection technique (see EalTUI 
above). In this work they show that the right inertia pa- 
rameter of the collective motion associated with the linear 
momentum operator P (N in our case) is obtained when 
the GCM coordinates are the position {ip in our case ) and 
the velocity (/i in our case). That means the dynamics as- 
sociated with Eq. 111! has the right inertial parameter. In 
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a VAP approach one always obtain the right mass param- 
eter in]. 

5.2 Pairing correlations. 

In a canonical ensemble the BCS order parameter is iden- 
tically zero. For this reason it is necessary to define an- 
other quantity to characterize pair correlations in a state 
of fixed numbers of electrons. We choose the pairing pa- 
rameter used in ref. ^Hl 

Z\b-G^Cfc (22) 
fc 

where the subindex b indicates the number parity of the 
grain. The quantities Cfe's are defined by 

Cfc = (4+Cfc+4-Cfc-) - {cl^Ck+){cl_Ck-) (23) 

and form a set of correlators which measure the fluctua- 
tions in the occupation numbers. The expectation values 
( ) are to be calculated with the wave functions of the 
corresponding approach using the formula developed in 
Appendices A and B. In an uncorrelated or in a blocked 
state one has Ck = 0. In the grand-canonical case the C^'s 
reduce to Ck — UkVk and Ab coincides with the usual su- 
perconducting order parameter. 

In Fig. [7|we show our results for the pairing parameter 
in units of A for even (upper panel) and odd (lower panel) 
systems respectively. As we can see in both plots the sharp 
transition occurring in the BCS and PBCS methods is ab- 
sent in the GCM approaches as well as in the exact solu- 
tion. The peculiar behaviors of Aj, in the exact and GCM 
approximations before and after the BCS breakdown are 
related to the change of a pairing delocalized in energy 
(weak pairing regime) to a localized one (strong pairing 
regime). The rough decrease of Ai, with N is connected to 
the special feature of the model, for which the constant G 
of Eq. (|22|l is inverse proportional to the number of elec- 
trons. The fact that A}, converges monotonically to the 
final value A , in the even case from above and in the 
odd one from below, is due to the blocking effect. In these 
plots we observe again that the GCMVAP approaches pro- 
vide solutions closer to the exact one than the GCMPAV 
approaches. 

5.3 Collective wave functions. 

We now look at the structure of the GCM states in the 
space of the collective parameter ^. The collective weights 
/(^) can not be interpreted as probability amplitudes be- 
cause the generating states \BCS{£,))n are not, in gen- 
eral, orthogonal to each other. The amplitudes G{£,) of 
Ea. H21() . on the other hand, play the role of "collective 
wave functions" , they are orthogonal and their modules 
squared have the meaning of a probability. 

The quantities |^/(0P are plotted in Figs. |H1 El and 
IIUI as a function of the parameters G, fi and AN"^, re- 
spectively. The behavior of |5(0P as a function of ^ indi- 
cates which are the most relevant components of the states 




80 160 240 320 400 480 560 
N 



Fig. 7. The gap parameter At for the different calculations as 
a function of the number of electrons for even and odd grains. 
For particle numbers smaller than those shown in the BCS plot 
the BCS gap parameter goes sharply to zero. 



\'^n{£.)) terms of the parameter To guide the eye 
we have also plotted in these figures the projected energy 
i?jv(f) of Eq. (|HJl. For simplicity we restrict our discussion 
to even systems. In Fig. |S1 we represent these quantities 
for the coordinate Gtriai in the GCMPAV and GCMVAP 
approaches. Since the wave functions of both approaches 
do not differ qualitatively we shall discuss both cases to- 
gether. The fact that the projected energies are lower in 
the PAV than in the VAP approach for a given Gtriai — Gc 
is due to the fact that Gc is almost zero for all particle 
numbers in the VAP approach while it varies considerably 
with the particle number for the PAV case, see Table |21 
We find broad potential energy curves for small particle 
numbers and narrower ones with increasing N. Interest- 
ingly the potential energy curves for the GCMPAV and 
GCMVAP approaches are rather different for small parti- 
cle numbers and they become similar for the large ones. 
Concerning the wave functions for = 20 we obtain very 
broad distributions corresponding to a situation of weak 
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Fig. 8. The projected energies En{G) versus G in the GCM- 
PAV (thin continuous hues) and GCMVAP (thick continuous 
Unes) approaches for even systems. The coUective wave func- 
tions \G{G)\'^ for the GCMPAV (thin dashed lines) and the 
GCMVAP (thick dashed lines) approaches in arbitrary units. 
The vertical scale applies for En{G), the minimum of En{G) 
in each approach has been set equal to zero. The top scale 
applies only for the top panel. 



pairing dominated by fluctuations in the order parameter 
A, see also panel a) of Fig. [3 For TV ~ 40 we find that the 
wave functions are not that extended anymore but they 
still present a two peak distribution, wit the first peak 
around the non-superconducting solution and the other 
around a superconducting one. For larger particle num- 
bers {N ^ 86,172,400) a one peak distribution emerges 
with the width of the peak getting smaller for increasing 
particle number. Looking at Fig. |21i we see that at large 
TV the distribution peaks around the wave function with 
a gap very close to the bulk one. 

The results for the /j. coordinate are presented in Fig.|5| 
For TV = 20, in the weak pairing region, we obtain a flat 
potential which shape corresponds to the physics already 
discussed in relation with Fig. O The collective wave func- 
tion, also according to the discussion of Fig.Eb and Fig. 13 
displays an oscillating behavior with maxima around fi = 
kd and minima (zeroes) around /i = k{d+l/2). The height 
of the maxima decreases considerably for k values different 
from 10 and 11, i.e., the collective wave function is mainly 
formed by the HF solution around the Fermi level and the 



-2 -1 1 -0.5 0.0 0.5 1.0 

^ - ^Bcs (meV) n - n'^cs (meV) 

Fig. 9. The projected energies i5jv(/i) (continuous lines) and 
the collective wave functions |C/(a*)|^ ( dashed lines) versus fi. 
The vertical scale applies for En{i-i-), the minimum of i5jv(/i) 
has been set equal to zero. |0(/i)P is in arbitrary units. The 
top scale applies only for the top panel. 



N=2G components of the BCS solution of the levels above 
and below the Fermi level. For TV = 40 the weak pairing 
regime persists and the wave function displays a structure 
similar to the TV = 20 case but with the strength much 
more concentrated owing to the fact that the level spacing 
decreases with increasing number of particles. For TV = 86 
the two peak structure is just a reminiscence of the weak 
pairing situation and for TV = 172 and 400 a one peak 
structure emerges indicating the strong pairing situation. 
The potential energy curves get steeper with increasing N 
and the localization of the peak around n^^g sharpens in 
the same way. As one can see in Fig.l^J) the range of the A 
parameter covered by the wave functions diminishes with 
increasing TV. 

Lastly we discuss the Z\TV^ coordinate in Fig. ^| The 
potential energy curves are easy to understand. In the 
small particle number limit the BCS solution does not pro- 
vide a superconducting solution and therefore (Z\TV^) = 0. 
On the other hand the BCS approximation provides the 
exact solution in the bulk limit, i.e. there (Z\TV^) >> 1. 
Accordingly we expect minima in the projected energies 
at small Z\TV^ for low N and at large AN"^ for large N. The 
potential energy curves get softer with growing N because 
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Fig. 10. The same as in figure|H|but for the parameter AN'^. 



and bulk regimes. Concerning the GCMPAV calculations 
we find that all three proposed coordinates, in spite of not 
being able to reproduce the exact results, describe qualita- 
tively the correct physics washing out the phase transition 
found in the BCS and the PBCS approaches. Concerning 
the degree of accuracy we find that the n coordinate is the 
most effective of all the three followed by the AN"^ one. 

We think that these results are rather general and ap- 
ply to many more complex Hamiltonians than the naive 
pairing one considered here. Since the GCM Ansatz in- 
cludes explicitly fluctuations in the wave function it is 
very appropriate to deal with finite systems where phase 
transitions may take place. The method, contrary to other 
approximations, applies equally well to systems with very 
few or very large particle number. On the other hand the 
GCM Ansatz is very versatile to be adapted to other phys- 
ical situations by considering additional coordinates to the 
ones discussed in this work. 

This work has been supported in part by DGI, Ministe- 
rio de Ciencia y Tecnologfa, Spain, under Project FIS2004- 
06697. M.A.F. acknowledges a scholarship of the Programa 
de Formacion del Profesorado Universitario (Ref. AP2002- 
0015) 



A Peculiarities of the PBCS approximation 

In this appendix we discuss some numerical aspects of the 
solution of the equations used in this article. 



for increasing AN"^ it is energetically easier to change this 
value. As it should be the potential energy curve for the 
GCMVAP approach lies below the GCMPAV one. Con- 
cerning the collective wave functions, their behaviors cor- 
respond to the shape of the potentials. For < 86 there 
is a finite probability of having an uncorrelated HF solu- 
tion as a component of the collective wave functions and 
only for A^ > 172 we obtain Wigner-like functions with 
zero amplitude for the HF component. The range of Z\'s 
covered by the wave functions can be read from Fig. 13;. 

Finally, we would like to mention that a very detailed 
comparison of our wave functions and the exact ones has 
been made in ref. ^Sl- We find that the physical content 
of the GCMPAV wave functions and the exact ones is 
identical. 



A.l Odd particle number case 

The PBCS and HW methods can be extended to systems 
with an odd number of particles by blocking one of the 
available states. A system with an odd particle number is 
described by the state 

X -[[(e-'f/^, + e'^^'vkcl^cljn (24) 

with N an even number and I the blocked state. The 
residuum integral of Eq. Q now looks like 



6 Conclusions 

In this paper we have presented a detailed formulation 
of the particle number projected Generator Coordinate 
Method. We have discussed two different coordinates to 
generate wave functions for the variation after projection 
method and three for the projection after variation one. 
The theory has been applied to study superconducting 
grains with a pairing Hamiltonian. We have shown that 
the GCMVAP calculations with both proposed coordi- 
nates reproduce the exact results in the weak, crossover 



J A/ 



1 

2^ 



dip e 



-i{M-2m)ip/2 



n 

,jM,l 



(25) 



in an obvious notation. As before all expectation values 
can easily be calculated in terms of the residuum integrals, 
for example the norm matrix element is given by 



N+i{BCS\BCS) 



N + l 



= Rn 



(26) 
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and the projected energy by 



E' 



7V+1 



N 



G 



2 } 



G 



N 

ujVjUkVk' 



The superindex is somewhat superfluous and it can 
be suppressed. Blocking different states I one obtains dif- 
ferent excited states. The lowest of these energies corre- 
sponds to the ground state of the system with an odd 
particle number. The PBCS and HW variational equa- 
tions obtained with theses states can be guessed without 
further calculations: In all sums and products in the equa- 
tions presented for even systems the term corresponding 
to the blocked state I has to be excluded. In the same way 
one can write down states with 3,5... etc. blocked states. 



the numerical solution of H30|l in an efficient way requires 
the computation of the residuum integrals in a fast and 
accurate way. For this purpose we have implemented Fast 
Fourier Transform routines to evaluate the integrals re- 
ducing the number of different integrals as much as possi- 
ble to minimize the computational effort. For this purpose 
the following two identities can be used. The first one was 
found by Dietrich et al jlTj. It can be shown that the 
residuum integrals satisfy the following recursion relation 



2 T>jl,---,jM,k 



vtR: 



(31) 



The knowledge of two residuum integrals allows to calcu- 
late a different one. This reduces the number of numerical 
integrations by one third. A second, more powerful rela- 
tion was found by Ma and Rasmussen j2()| . 



1 



2(A/-m-l) 2m 



J A/ 



J=3l, 



■ jAf 



A. 2 Numerical solution of the PBCS equations. 

The set of coupled non- linear equations Q is usually solved 
by an iterative procedure. In order to speed up such proce- 
dure it is convenient to introduce the variable Xk through 
the relation 



Xk 



1 



1 



Xk 



1 + 



1 has been 



where the normalization condition uf. + v 
taken into account by construction. In terms of Xk the 
variational equations look like 



1 /2 

{ik + Ak)Xk ~^k{xk 
The additional transformation 

Xk = exp dk 



1) 



0. 



(28) 



(29) 



allows to isolate the new variable in terms of the fields 
ik^Ak and Ak 



n 



n.Q. 



(32) 



This formula allows to calculate all residuum integrals 
if the integrals i?Q and i?Q are known. This relation reduces 
to + 1 the overall number of numerical integrations for 



(27) a given set of Vj's and Mj's. 



In the PBCS one only needs Ma's relation to calculate 



three terms: 



- R{^^ and R. 



.jk 



i?f . If R° and 



Rq are known, i?j can be obtained by Dietrich's recursion 
relation. First we consider i?^ with m = 1 or m = 2. Ma 
and Rasmussen's formula reduces to 



+ i-y 



/-m,,2 tdJ /-m,,2 jDk 



(33) 



where we have used the identity Qj 



1 



1 . The difference 



i?f - Ri" can be written in a simplified way as follows 



jjk 



= 2 sinh" 



A, 



2Ak 



(30) 



The right hand side of this equation does not depend 
explicitly on Ok because the fields ek,Ak and Ak are in- 
dependent of Ok- This fact is very useful to solve Eq. 
by numerical iteration. We start with a guess of 9k (for 
instance the grand-canonical solution) and solve un- 
til the convergence of the energy, to a given tolerance, has 
been reached. 

The variational equations involve the computa- 
tion of the residuum integrals R'^^'"'^" . These integrals 
can be calculated analytically using the existing closed 
analytical expression j25) . Their evaluation, however, re- 
quires the addition of many terms making the whole com- 
putation a very time consuming approach. As the integrals 
must be computed many times for different sets of Uk,Vk, 



R 



.jk 



R{ 



jk 



(34) 



The calculation of i?^*^' — R{'^'' is a bit more compli- 
cated, the result is 



j^jki _ ^ki ^ 



vf) 



{vl - v'^){vl - vf) 



{vf -ij){vf -vD' 



(35) 



Since the indices of the residuum integrals can be per- 
muted, i?^*^' — R{''^ can be expressed for all possible com- 
binations of Wj, Ufe, vi by the equations above. 
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